A high-quality chromosome-level genome assembly of the Chinese medaka Oryzias sinensis

Oryzias sinensis, also known as Chinese medaka or Chinese ricefish, is a commonly used animal model for aquatic environmental assessment in the wild as well as gene function validation or toxicology research in the lab. Here, a high-quality chromosome-level genome assembly of O. sinensis was generated using single-tube long fragment read (stLFR) reads, Nanopore long-reads, and Hi-C sequencing data. The genome is 796.58 Mb, and a total of 712.17 Mb of the assembled sequences were anchored to 23 pseudo-chromosomes. A final set of 22,461 genes were annotated, with 98.67% being functionally annotated. The Benchmarking Universal Single-Copy Orthologs (BUSCO) benchmark of genome assembly and gene annotation reached 95.1% (93.3% single-copy) and 94.6% (91.7% single-copy), respectively. Furthermore, we also use ATAC-seq to uncover chromosome transposase-accessibility as well as related genome area function enrichment for Oryzias sinensis. This study offers a new improved foundation for future genomics research in Chinese medaka.


Background & Summary
Chinese medaka, or Oryzias sinensis, is a teleost fish closely related to Japanese medaka (Oryzias latipes), which has been used as a model organism in many genomic studies as well as in aquatic toxicology research 1 .Both fish belong to the family Adrianichthyidae, commonly referred to as the Medaka family.Similar to its relative, Chinese medaka is also attracting the attention of scientists due to its small size and short generation interval 2 .There are several differences between these two fishes, including their vertebrae and pectoral fin strip; most importantly, Chinese medaka is mainly distributed in freshwater while Japanese medaka can adapt to a certain level of salinity.Because of this, Chinese medaka could be more suitable for freshwater quality evaluation.
Although the Chinese medaka genome has been released 3 , the completeness and genome annotations still need to be further improved.The reported genome was only released at the scaffold level with alignments to the chromosomes of Japanese medaka.Several phylogenetic analyses and taxonomic revisions for medaka have already found that the Chinese medaka is different from its Japanese relative in chromosome constitution; the diploid chromosome number is 46 in O. sinensis, and 48 in O. latipes 4,5 .Therefore, a high-quality reference genome for O. sinensis is increasingly important to support future research.
In the present study, an improved high-quality chromosome-level genome assembly of O. sinensis was generated using single-tube long fragment read (stLFR) reads, Nanopore long-reads, and the Hi-C sequencing data.The genome size was 796.58 Mb with a scaffold N50 length of 30.38 Mb (Table 1).A total of 712.17 Mb (89.34%) of assembled sequences were anchored to 23 pseudo-chromosomes, with lengths ranging from 58.48 Mb to 18.87 Mb (Table 2).Based on this improved genome assembly, repeat elements and gene structure annotations were conducted combining de novo prediction, homolog-based alignment, and transcriptome-assisted methods.
Benchmarking Universal Single-Copy Orthologs (BUSCO) evaluation result showed that the final assembly was benchmarked at 95.1% and the annotation reached 94.6% (Table 3).
In addition, we conducted ATAC-seq try to find out the chromosome accessibility of O. sinensis.After we got the peak results we also performed function enrichment of related genome areas and obtain some clue of transcriptional activity.Taken together, this study could provide a new reliable foundation for research on the Chinese medaka, as well as for its use as a model organism.

Methods
De novo genome assembly.The samples of O. sinensis were obtained from the National Plateau Wetland Research Center, Southwest Forestry University; muscle tissue was used for nucleic acid extraction.High-quality purified RNA was used to construct a transcript sequencing library and DNA was used to construct stLFR, Nanopore long-reads, Hi-C sequencing, and ATAC-seq libraries.stLFR technology added the same barcode sequence to subfragments of the original long DNA molecule, allowing these co-barcoded sub-fragments to be subsequently sequenced using a second-generation platform 6 .This method generates long reads with high accuracy, enabling high-quality assembly and subsequent analysis.Most importantly, stLFR was cost-effective and has been widely used in aquatic genome projects 7,8 .The genome size and heterozygosity of O. sinensis were estimated using Jellyfish v2.2.6 9 and GenomeScope v1.0.0 10 by k-mer analysis with clean stLFR data.The clean data was then used in de novo genome assembly with the stlfr2supernova pipeline (https://github.com/BGI-Qingdao/stlfr2supernova_pipeline).This pipeline conducts de novo assembly using stLFR data with Supernova Assembler, which refers to the de novo software from 10X Genomics 11 .Nanopore long-reads were employed to carry out further scaffolding, gap-closing, and polishing at the same time using TGS-GapCloser 12 .The size of the assembly after these steps was bigger than the k-mer result (approximately 906 Mb) so purge_haplotigs 13 was applied to remove redundancy.Finally, a draft assembly that covered approximately 796 Mb of the genome with a contig N50 length of 373.71 Kb was obtained.
analysis and chromosome assembly.The Hi-C reads were aligned to the draft assembly using HiCPro 14 to find valid read pairs.Juicer 15 and 3D-DNA 16 were then used to finish the construction of chromosomes, with manual correction of misjoins, wrong order, and opposite orientation using Juicebox 15 .The Hi-C scaffolding resulted in 23 pseudo-chromosomes with a total length of 712.17 Mb.We offer a chromosome circus map using TBtools 17 and heatmap using Juicebox 18 (Fig. 1A,B).Genome assembly statistics are shown in Table 1 and Table 2. BUSCO (Benchmarking Universal Single-Copy Orthologs) 19 evaluation of the assembly reached 95.1% using actinopterygii_odb9 database, indicating a well-covered fish genome assembly (Table 3).

Repeat annotation.
Tandem repeats and interspersed repeats identification was essential before protein-coding gene and function annotations.RepeatMasker 20 and RepeatProteinMask were applied to find repeat elements as homology predictions based on RepBase 21 .For the de-novo method, RepeatModeler 22 was used to predict repeat elements; LTR_FINDER 23 , and TRF tool 24 were also used to predict repeat elements based on sequences features.Overall, 311.32 Mb of the O. sinensis genome assembly were identified as repetitive elements, accounting for 39.05% of the whole genome (Table 4).
Gene prediction and annotation.Protein coding genes were predicted with multiple sources of evidence including homology-based alignments, de novo prediction, and transcriptome-assisted methods.In homology alignments, the protein sequences of Danio rerio (GRCz11), Ictalurus punctatus (IpCoco_1.2),Oryzias javanicus (OJAV_1.1),Oryzias latipes (ASM223467v1), Oryzias melastigma (Om_v0.7),Poecilia formosa (PoeFor_5.1.2),and Takifugu rubripes (FUGU5) were mapped to a repeat soft masked draft genome using blat 25 , and Genewise was applied to define gene models 26 .In de-novo prediction, Augustus 27 was used to predict the coding regions of genes.In transcriptome-assisted methods, two different methods were used to obtain predicted gene sets.First, RNA reads were mapped to the genome assembly using Hi-SAT2 28 ; the result was used to build a transcript gene model using Stringtie 29 and TransDecoder.Second, a transcriptome was first assembled with RNA reads using Fig. 1 (A) Global genome landscape of Chinese medaka.From inner to outer circles: Density of genes with 500 kbp windows, ranging from 0 to 50; Density of TE with 500 kbp windows, ranging from 0 to 600; Density of TRF with 500 kbp windows, ranging from 0 to 450; Distributions of GC content with 500 kbp windows.(B) Hi-C interaction heat map and scatter plot of genome assembly GC content and sequencing depth.
Trinity, before PASA was used to get a gene model.Finally, all evidence was merged to form a consensus gene structure annotation result using GLEAN 30 .A total of 22,461 protein-coding genes were successfully predicted.BUSCO 19 assessments reached 94.6% using the actinopterygii_odb9 database, indicating relatively complete gene annotation coverage (Table 3).
All predicted genes were compared with public biological function databases including KEGG 31 , Swissprot 32 , and TrEMBL 33 (https://www.uniprot.org/statistics/TrEMBL)using BLASTp 34 with an E-value cutoff of 1e-5 for functional annotation.Interpro was applied to provide functional analysis of protein sequences by classifying them into families and predicting the presence of domains and important sites, followed by Gene Ontology annotation [35][36][37] .Overall, 22,162 protein-coding genes (98.67%) were successfully functionally annotated (Table 5).
Non-coding RNAs (ncRNAs) are an important part of genome annotation as they can be active in transcriptional and translational regulation of gene expression as well as in the modulation of protein function 38 .Since ribosomal RNA (rRNA) was highly conservative, vertebrate rRNA data was used as a reference to map to the draft genome using BLASTn with an E-value of 1e-5.For the discovery of transfer RNA (tRNA), tRNAscan-SE v1.3.1 was applied with eukaryotic parameters according to the characteristics of tRNA 39 .Rfam database was applied to find microRNAs (miRNA) and small nuclear RNA (snRNA) 40,41 .In total, 676 rRNA, 740 tRNA, 231 miRNA, and 1,168 snRNA genes were identified from the O. sinensis genome.
Chromosome accessibility analysis using ataC-seq.ATAC-seq is a widely used method used to identify regions of open chromatin in genomes.SOAPnuke (v1.5.2) 42 was used to filter out reads with low quality original data or high content of unknown bases.Bowtie2 (v2.2.5) 43 was used to align reads to the O. sinensis genome assembly and peak calling was performed with MACS2 software (v2.1.0) 44.Peaks related to genes were used to carry out GO and KEGG enrichment analysis (Fig. 2).

Data Records
The Nanopore long reads, stLFR genomic sequencing data, Hi-C data as well as the assembled genome have been deposited at China National GeneBank DataBase (CNGBdb) under the accession CNP0003475 45.The whole sequencing dataset of O. sinensis was deposited in the NCBI Sequence Read Archive database (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA895195)under project identification number PRJNA895195.This Whole Genome Shotgun project has been deposited at DDBJ/ENA/GenBank under the accession JAUDJI000000000.The version described in this paper is version JAUDJI010000000 46 .Sequence Read Archive (SRA) project number was SRP410304 47 .DNA sequencing data from the WGS library were deposited in the SRA at SRR22435960 48 .DNA sequencing data from the Hi-C library were deposited in the SRA at SRR22435961 49 .The ATAC-sequencing data were also deposited at Figshare 50.

technical Validation
To evaluate the quality of the genome assembly, stLFR reads were mapped to the final reference genome assembly using BWA (v0.7.12).A total of 98.27% of reads were mapped, covering 98.29% of the genome sequence.Genome average sequencing depth reached 172x, and 97.34% of the genome had a sequencing depth of over 20x.The scatter plot between assembly sequencing depth and GC content found no abnormal GC content, recognized as no exogenous pollution.The completeness of the genome assembly and annotation was assessed using BUSCO (v3.0) with the actinopterygii_odb9 database.The BUSCO benchmark of genome assembly and gene annotation reached 95.1% and 94.6%, respectively.

Table 1 .
Genome assembly statistics of Chinese medaka Oryzias sinensis.

Table 2 .
Chromosome length distribution of assembled genome.*Unplaced referring the sequences which did not mount to known chromosome.

Table 3 .
Genome and Annotation BUSCO evaluation result.

Table 4 .
Prediction of repeat elements.

Table 5 .
Genome function annotation result.